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ABSTRACT 

DTIC  QUALITY  INBPECTED  9 

We  consider  preconditioning  methods  to  accelerate  convergence  to  a  steady  state  for  the 
incompressible  fluid  dynamic  equations.  The  analysis  relies  on  the  inviscid  equations.  The 
preconditioning  consists  of  a  matrix  multiplying  the  time  derivatives.  Thus  the  steady  state 
of  the  preconditioned  system  is  the  same  as  the  steady  state  of  the  original  system.  We 
compare  our  method  to  other  types  of  pseudo-compressibility.  For  finite  difference  methods 
preconditioning  can  change  and  improve  the  steady  state  solutions.  An  application  to  viscous 
flow  around  a  cascade  with  a  non-periodic  mesh  is  presented. 
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1  Introduction 


One  way  to  solve  the  steady  state  incompressible  equations  is  to  march  the  time  dependent  equations 
until  a  steady  state  is  reached.  Since  the  transient  is  not  of  any  interest  one  can  use  accderation 
techniques  which  destroy  the  time  accuracy  but  enables  one  to  reach  the  steady  state  faster.  Such 
methods  can  be  considered  as  preconditionings  to  accelerate  the  convergence  to  a  steady  state. 
For  the  incompresible  equations  the  continuity  equation  does  not  contain  any  time  derivatives. 
To  overcome  this  difficulty,  Chorin  [3]  added  an  artificial  time  dmivative  of  the  pressure  to  the 
continuity  equation  together  with  a  multiplicative  variable,  $  .  With  this  artificial  term  the  resultant 
scheme  is  a  symmetric  hyperbolic  system  for  the  inviscid  terms.  Thus,  the  system  is  well  posed 
and  and  numerical  method  for  hyperbolic  systems  can  be  used  to  advance  this  system  in  time.. 
The  free  parameter  is  then  chosen  to  reach  the  steaidy  state  quickly.  Later  Turkel  ([8],  [9],  [10]) 
extended  this  concept  by  adding  a  pressure  time  derivative  also  to  the  momentum  equations.  The 
resulting  system  after  preconditioning  is  no  longer  symmetric  but  can  be  symmetrized  by  a  change 
of  variables. 

Thus,  we  wUl  consider  systems  of  the  form 

to*  +  /*  +gy  =  0. 

This  system  is  written  in  conservation  form  though  for  some  applications  this  is  not  necessary.  Our 
analysis  will  be  based  on  the  linearized  equations  so  the  conservation  form  does  not  appear  in  the 
analysis  though  it  does  appear  in  the  final  numerical  approximation.  This  system  is  now  replaced 
by 


P"*to«  + /,  +  py  =  0,  (1) 

or  m  Unearized  form 

+  AWf  +  Bwy  =  0,  (2) 

with  A  and  B  constant  matrices. 

For  this  system  to  be  equivalent  to  the  original  system,  in  the  steady  state,  we  demand  that 
P~^  have  an  inverse.  This  only  need  be  true  in  the  flow  regime  under  consideration.  We  shall 
see  later  that  frequently  P  is  singular  at  stagnation  points.  Thus,  we  will  temporarily  consider 
strictly  flows  without  a  stagnation  point.  We  also  assume  that  the  Jacobian  matrices  A  = 
and  B  =  are  simultaneously  symmetrizable.  In  terms  of  the  ‘symmetrizing’  variables  we  also 
demand  that  P  be  positive  definite.  We  shall  show  later,  in  detail,  that  it  does  not  matter  which 
set  of  dependent  variables  are  used  to  develop  the  preconditioner.  One  can  transform  between  any 
two  sets  of  variables.  Thus,  when  we  are  finished  we  will  analyze  a  system  which  is  similar  to  (2), 
where  the  matrices  A  and  B  are  symmetric  and  P  is  both  symmetric  and  positive  definite.  Such 
systems  are  known  as  symmetric  hyperbolic  systems.  One  can  then  multiply  this  system  by  w  and 
integrate  by  parts  to  get  estimates  for  the  integral  of  wif,  i.e.  energy  estimates.  These  estimates 
can  then  be  used  to  show  that  the  system  is  well  posed  .  We  stress  that  if  P  is  not  positive  then 
we  may  change  the  physics  of  the  problem.  For  example,  if  P  =  —  /  then  we  have  reversed  the  time 
direction  and  must  therefore  change  all  the  boundary  conditions.  Keeping  the  right  signs  for  the 
eigenvalues  is  a  necessary  but  not  sufficient  condition  for  well-posedness. 

With  this  assumption  the  steady  state  solutions  of  the  two  systems  are  the  same.  Assuming 
that  the  steady  state  has  a  unique  solution,  it  does  not  matter  which  system  we  march  to  a  steady 
state.  We  shall  later  see  that  for  the  finite  difference  approximations  the  steady  state  solutions  are 
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not  necessarily  the  same  and  nsnally  the  preconditioned  system  leads  to  a  better  behaved  steady 
state. 

2  Incompressiblle  Equations 

Consider  the  incompressible  inviscid  equations  in  primitive  variables. 

ttx  +  Uy  =  0 

+  Uttx  +  VUy  +px  =  0 

+  VVy  +Py  =  0 

We  generalize  Chorin’s  pseudo-compressibility  method  [3].  Using  the  preconditioning  suggested  in 
[8]  (with  a  =  1)  we  have 

1 

■^Pt  +  «x  +  t>y  =  0 

U 

-^Pt  +  tit  +  ttttx  -b  t>«y  -i-  Px  =  0 

V 

-^Pt  +Vt  +  UVt  +  VVy  +  Py  =  0  (3) 

or  in  conservation  form 

0 
0 
0 


D  =  u;%A -h  ufjB  —  1  < <  1 


^P,  +  Ur  +  V,  = 

+  »!  +  («’ +?)»  +  («<>),  = 
2v 

■^Pt  +  Vt  +  («w),  +  (u*  +  p),  = 

We  can  also  write  (3)  in  matrix  form  using 

/  0  0  \ 
p-»  =  u/0^  10, 

\v/0^  0  1/ 

/  0  0  \ 

P  =  -ti  1  0 

\  -r  0  1  / 

Multiplying  by  P  we  rewrite  this  as 

wt  +  Pi4ti>x  +  "PBwy  =  0. 
Let 
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where  idfi,u)2  are  the  Fourier  transform  variables  in  the  x  and  y  directions  respectively.  The  speeds 
of  the  waves  are  now  governed  by  the  roots  of  det(XI  -  PAuJi  —  PBuf^)  =  0  or  equivalently 
det{XP~^  —  Auii  —  BcJj)  =  0.  Let 


q  =  tia»i  +  viJ2- 


Then  the  eigenvalues  of  PD  are 


do  =  g  (4) 

d±  =  ±/J 

and  so  the  ‘acoustic*  speed  is  isotropic. 

The  spatial  derivatives  involve  symmetric  matrices,  i.e.  D  is  a  symmetric  matrix  but  P  is  not 
symmetric.  Thus,  while  the  original  system  was  symmetric  hyperbolic  the  preconditioned  system 
is  no  longer  symmetric.  In  [8]  it  is  shown  that  as  long  as 

/3^  >  (tt*  +  »*) 

the  equations  can  be  symmetrized.  On  the  other  hand  the  eigenvalues  are  most  equalized  if  = 
(it^  + [8].  So,  we  wish  to  choose  slightly  larger  than  u^  +  v^.  However,  numerous  calculations 
verify,  that  in  general,  a  constant  ^  is  the  best  for  the  convergence  rate.  The  reasons  for  this  are 
not  clear. 

We  wish  to  stress  that  has  the  dimensions  of  a  speed.  Therefore,  cannot  be  a  universal 
constant.  There  are  papers  that  claim  that  =  1  or  —  2.5  are  optimal.  Such  claims  cannot  be 
true  iir  general.  It  is  simple  to  see  that  if  one  nondimensionalizes  the  equation  then  gets  divided 
by  a  reference  velocity.  Hence,  the  optimal  ‘constant’  depends  on  the  dimensionalization  of  the 
problem  and  in  particular  depends  on  the  inflow  conditions.  In  many  calculations  the  inflow  mass 
flux  is  equal  to  1  or  alternatively  p  4-  (tt^  +  v^)/2  —  1.  Such  conditions  wiU  give  an  optimal  close 
to  one. 

We  next  define  the  Bernoulli  function 

ir  =  p  +  (t»»  +  r»)/2. 

Bernoulli’s  theorem  states  that  when  the  flow  is  steady  and  invisdd  then  H  is  constant  along 
streamlines.  We  now  multiply  the  second  equation  of  (3)  by  u  and  the  third  equation  of  (3)  by  v 
and  add  these  two  equations.  If  4-  the  result  is 

Ht  4-  4-  vH^  =  0.  (5) 

Thus,  by  altering  the  time  dependence  of  the  equations  we  have  constructed  a  new  equation  in 
which  H  is  convected  along  streamflnes.  FVirthermore,  if  H  is  a  uniform  constant  both  initially  and 
at  inflow  then  H  wUl  remain  constant  for  afl  time.  On  the  numerical  level  this  will  usually  not  be 
true  because  of  the  introduction  of  an  artificial  viscosity  or  because  of  upwinding.  For  viscous  fiow, 
(5)  is  replaced  by 


Ht  4-  uHx  4-  vHy  =  — (uAu  4-  uAu) 
H/t 
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We  note  that  these  relationships  for  H  follow  from  the  momentum  equations  and  do  not  depend 
on  the  form  of  the  continuity  equation.  Hence,  we  consider  the  following  generalization  of  (3) 


Pt  +  aHt  +  tig +  Vy  =  0 


ou 

av 


Pt  +  tlt+  UUg  +  VUy  +Pg  =  0 


•^Pt  +  Vt-^UVg  +  VVy+Py  =  0 

where,  a,a  and  /5  are  free  parameters.  When  u;]  +  <t;|  =  1  the  eigenvalues  of  PD  are 

s±y/?TWd 


(6) 


(Jltl  +  U>2Vt 


2d 


where 


d  =  1  +  a  -  o^,  s  =  (1  -  Q)(ufiu  +  Wju). 


Hence,  the  ‘acoustic’  eigenvalue  is  isotropic  if  a  =  1.  Furthermore,  d  =  1  if  either  a  =  0  or 
+  v^.  For  a  =  0  we  recover  our  ori^al  scheme.  For  o  =  — 1  the  time  derivative  of  the 
pressure  no  longer  appears  in  the  continuity  equation.  For  general  a,  /3  we  have 


,  /  (o+ 1)  au  av  \ 
^  au  0  , 

^  \  av  0  0^} 


P  = 


/  — ou 

-au  1  +  0- 
—av 


1  +  a  —  oa 


If  we  write  the  equation  in  conservation  form  (1)  we  have 


P“* 

eoiMervaUve 


1  /  (o+l 

^  V(l  +  a  + 


—av 

aotnv 

“gS” 

1+0- 


1)  ou 

a)u  0^  +  ou*  out? 
a)u  ouv  0^ 


av  \ 
auu  I  , 

+  OU*  J 


coruervative  — 


1  +  0- 00*2^  * 


^  +  o(u*  +  u*)  -ou 

-(1  +  o  +  a)u  1  +  o  -  2^ 
^  -(l  +  o  +  a)u  2^ 


— ou 

aoniv 

1+0-22^ 


In  [9]  an  analogy  to  the  symmetric  preconditioning  of  van  Leer,  Lee  and  Roe  was  constructed 
for  the  incompressible  equations,  ff  we  choose  o  =  1  P  is  symmetric.  If  we  also  choose  ,3^  =  u*  +  u* 
then  we  get  the  preconditioning  of  van  Leer  et.al.  . 


u*+u* 


P  = 


-u  l  +  -st 


— u 


—V 

uv 


— U 

a 

yv  _  1 1  .y  „ 
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These  examples  show  that  preconditioning  is  not  unique.  If  fact,  since  the  determinant  of  the 
transpose  of  a  matrix  is  equal  to  the  determinant  of  the  ori^al  matrix  it  follows  that  the  transpose 
of  P  is  also  a  preconditioner  with  the  same  eigenvalues  for  the  preconditioned  system.  These  various 
systems  will  have  the  same  mgenvalues  but  different  eigenvectors  for  the  preconditioned  system. 
Numerous  calculations  show  that  the  system  given  by  P  in  (3)  is  more  robust  and  converges  faster 
than  with  the  transpose  preconditioner.  This  shows  that  it  is  not  sufficient  to  consider  just  the 
eigenvalues  but  that  the  mgenvectors  are  also  of  importance.  The  eigenvectors  are  given  in  ([10]). 
We  next  consider  the  preconditioner  considered  by  de  Jonette,  Viviand  et.  al.  ([4]).  Define: 

q  =  utt, +  t>tt,+p*-(r^  +  T*^) 
r  =  uvg  +  vv^+jt^-  (r^  +  t„) 
s  —  uq  +  vr 

Then  they  consider  the  following  extension  of  the  incompressible  Navier-Stokes  equations. 


Pt  +  dv(«s  +  Vy)  +  av»  =  6 

tt«  +  avq  +  +  «»)  +  evv*  =  0  (7) 

Vi  +  ayr  +  +  Vy)  +  ev»s  =  0 

In  the  steady  state  q  =:  r  =:  s  =  0  and  ««  +  Vy  =  0  and  so  we  recover  the  usual  incompressible 
equations.  av,dv,evtOtVffiv  are  free  parameters  that  satisfy  the  frdlowing  conditions 


av^  ~  dyev 
ov  >  0  dy  >9 
{dy  +  ayU^){dy  +  /Jyl/*)  >  0 


In  addition,  in  order  for  the  speed  of  the  convective  wave  to  remain  unchanged  we  add  the  condition 
ay  =  1.  FYom  the  momentum  equations  we  obtain 

[uu«  -hwi  +  ffyU^jUg  +  Vy)] 
l  +  eyU^ 


H«ice  we  can  rewrite  (8)  as 


^  ^  Pt  - 

+11,  +  Vy  =  0 

/8v«  _ 


/JyV 


p«  +  V,  +  r 


Comparing  this  with  (fi)  we  see  that  the  two  approaches  are  identical  if 
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dv 

ay 


«v 


(o+l)^-oat;2 


aa 


dy 


Choosing  a  =  1  and  0  =  we  get  the  standard  preconditioning  (3).  The  Viviand  parameters 
become  ay  =  — o,  fiy  —  —1,  ay  =  1,  dv  =  ey  =  a  .  Then  o  =  0  gives  the  Torkel  preconditioner 
and  a  =  1  ^ves  the  van-Leer  («ymmetric)  preconditioner. 


3  Difference  equations 


Until  now  the  entire  analysis  has  been  based  on  the  partial  diiferential  equation.  We  now  make 
some  remarks  on  important  points  for  any  numerical  approximation  of  this  system.  When  using 
a  scheme  based  on  a  Riemann  solver  this  solver  should  be  for  the  preconditioned  system  and  not 
the  original  scheme.  When  using  a  central  difference  schemes  there  is  a  need  to  add  an  artificial 
viscosity.  Accuracy  is  improved  for  low  Mach  number  flows  if  the  preconditioner  is  applied  only 
to  the  physical  convective  and  viscous  terms  but  not  to  the  artificial  viscosity.  The  use  of  a 
matrix  artificial  dissipation  ([7])  should  be  based  on  the  preconditioned  equations  as  for  Riemann 
solvers,  difference  scheme.  Hence,  both  for  upwind  and  central  difference  schemes  the  Riemann 
solver  or  artificial  viscosity  should  be  based  on  P~*|Pi4|  and  not  |i4|  i.e.  in  one  dimension  solve 
wt  +  P/t  =  P(P~^|PA|to«)«  .  When  using  characteristics  for  extrapolation  at  the  boundaries 
it  ?hould  be  based  on  the  characteristics  of  the  modified  system  and  not  the  physical  system. 
Preconditioning  is  even  more  important  when  uring  multigrid  than  with  an  explicit  scheme.  With 
the  orif^al  system,  the  stiffness  of  the  eigenvalue  greatly  affects  the  smoothing  rate  of  the  slow 
components  and  so  slows  down  the  multigrid  method,  [6].  We  conclude  that  the  steady  state 
solution  of  the  preconditioned  system  may  be  different  from  that  of  the  physical  system.  Thus,  on 
the  finite  difference  level  the  preconditioning  can  improve  the  accuracy  as  well  as  the  convergence 
rate. 

We  next  consider  adding  artificial  viscosity  to  the  system  (3).  We  first  rewrite  this  system 
eliminating  pt  from  the  velodty  equations.  This  ^ves 


or  in  matrix  form 


Pt  +  +  Vy) 

th+Px  +  VUy-  UVy 
Vt  +  UVg  -  VUs  +  Py 


-  0 
=  0 
=  0 


0  0  /  p  \ 

0  o  -«  tt  =0 

0  0  1  MWv 


(8) 
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wt  +  PAwg  +  PBwy  =  0. 

We  next  consider  the  use  of  a  matrix  valued  viscosity.  Let  D  =  uJiA  +  UiB  with  w?  +  =  1. 

The  non-preconditioned  matrix  viscosity  is  given  by  |i4|  in.  the  x  direction  and  |fi|  in  the  y  direction. 
Then 

(2q^  q  qRS 

q  R^  (Aj>  +  Aa^)  R^  +  (Aa  -  A3)  S^\R\  RS  (Aj^  +  Ag^  -  (Aj  -  A3) \R\) 
qRS  RS  (Aj^  +  Aa^  -  (Aj  -  A3)  |ii|)  (Aj*  +  Aa^)  5*  +  (Aj  -  A3) R^  |iJ| 


,  R  +  VWT4  ,  R-y/Wn 

Aa  = - 2 - ’  - 2 - 

R  =  tla^l  +  uwj,  5  =  VW2  -  VU2 

For  the  preconditioned  artificial  viscosity  we  consider  instead  P“'|PA|  and  P“^|PB|  (see  [7]). 
We  consider  the  case  a  =  1  with  0  and  a  arbitrary.  Then 

/  V12  Vi3  \ 

P-1|PD|  =  Vj,  V22  V23  . 

V  ^3  ^3  / 


l  +  a.«SVX 


V  .f  “  I 


V  “  I  Y  I 


,,  '  ,  -3  V  RS(?‘-a^)„X 


— ? - 


0  0‘X‘v‘X 

*'33  —  “7=  H - 5 - r 

Vd  r 
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where 


Va3  = 


(fi^Sv+Ru  (q^-0^))  O3^Su+Rv(0^-aq^))  X 

/3i,3 


V3J  = 


(fi^Sv  +  Ru  (aq^-0^))  Su-i- Rvift^  -  ^))  X 

fiiqi 


ttcJi  -ftwa  j  tifa^a  -t>a;i 


=  1  +  a  -  =  It*  + 

P* 


X  = 


mvd-p 

dZ 


j9* 


By  inspectioii  the  matrix  is  syHimetric  when  «  =  1.  For  the  special  case  o  =  1  and  /?*  =  «*  + 1>* 
the  fomialas  simfdify  and  we  get 


Vn  = 


A  + 1 


V21  =  Via  =  R- 
9 


V31  =  Vi3  =  R- 
9 


*  f  + 


e*(il  -  1)  _  a*  +  m* 

«  “  f 


Va3»VM 


«e(<-  1) 

f 


For  the  equations  in  conservatkm  form  we  midtiply  the  continuity  equation  by  n  and  add  to 
the  X  velocity  eqnation.  We  also  mnlti]^y  the  contmtdty  equation  by  v  and  add  to  the  y  vriocity 
equation. 


8 


4  Computational  Results 

We  BOW  present  a  calculation  for  two  dimensional  flow  around  an  cascade  to  demonstrate  the 
previous  theory.  The  discretization  is  based  on  the  multistage  time  method  coupled  with  a  central 
difference  approximation  as  described  in  ([5],  [7]).  The  basic  scheme  is  accelerated  by  using  a  local 
time  step,  residual  smoothing  and  multigrid.  This  code  was  further  developed  to  consider  cascade 
configurations  in  which  the  grid  is  not  necessarily  continous  across  the  wake  ([1],  [2]).  We  compute 
the  flow  about  a  NACA0012  with  periodic  external  boundaries.  The  flow  is  turbulent  and  we  use 
a  Baldwin- Lomax  turbulence  model,  with  Rt  =  500, 000,  Pr  =  0.7,  Prt  =  0.9  At  inflow  the  angle 
of  attack  is  specified  as  well  as  the  Bernoulli  constant,  p  +  =  1.  The  mesh  is  192  x  32  and 

is  shown  in  figure  1.  We  use  a  four  stage  Runge-Kutta  method  as  a  smoother  for  a  full  multigrid 
iteration.  We  choose  o  =  0  and  =  max(K(u^  -f-  with  K  =  1.1, =  0.4,  see  (6).  In 

figure  2  we  plot  the  convergence  rate  for  different  values  of  a.  We  see  that  the  fastest  convergence 
occurs  when  a  =  1  followed  by  a  =  0  and  finally  ot  =  —  1.  We  also  considered  viscous  flow  about 
a  VKI  cascade  (figure  3).  In  this  case  the  convergence  of  all  the  methods  slowed  down,  a  =  1  was 
stUl  the  most  efficient  method  but  the  differences  were  less  dramatic  than  in  the  previous  case.  In 
other  cases  in  was  necessary  to  choose  0  almost  constant.  The  symmetric  preconditioner,  a  =  1 
was  more  robust  but  not  faster  than  a  =  0. 

5  Conclusions 

A  three  parameter  preconditioning  matrix  has  been  introduced  for  the  incompressible  inviscid 
equations.  This  is  equivalent  to  the  pseudo-  compressibility  methods  considered  by  de  Jouette  et. 
al.  When  a  =  1  the  ‘acoustic’  speeds  are  symmetric.  Furthermore,  one  can  choose  the  parameter 
a  so  that  the  preconditioning  matrix  is  symmetric.  For  the  inviscid  case  considered  computed  a 
considerable  increase  in  the  convergence  rate  was  achieved. 

In  addition  the  incompressible  equations  offer  a  theoretical  advantage  over  the  compressible 
equations  for  the  theoretical  study  of  preconditioning  methods.  This  is  because  of  the  simpler 
nature  of  the  equations  and  the  fact  that  the  ori^al  method  of  Chorin  is  already  symmetric. 
Nevertheless,  a  central  difference  scheme  coupled  with  a  Runge-Kutta  time  advancement  suffers 
from  lack  of  robustness.  In  particular  needs  to  be  bounded  away  from  zero  at  a  relatively 
high  level  for  many  of  the  cases.  Using  the  symmetric  preconditioner  a  =  a  =  1  yields  a  more 
robust  scheme  though  it  does  not  seem  to  converge  faster  than  the  nonsymmetric  preconditioner. 
Furthermore,  changes  of  the  physical  inflow  boundary  condition  can  greatly  affect  the  choice  of  the 
optimal  a  and  /?.  The  major  increases  in  the  convergence  rate  are  for  the  Euler  equations.  For  the 
Navier-Stokes  equations  it  is  necessary  to  reformulate  the  preconditioning  matrix  to  account  for 
the  viscous  effecrs. 
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1.  Mesh  for  turbulent  flow  around  a  NACA0012  cascade. 
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